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Abstract 

The level-set method of topology optimization is used to design isotropic two- 
phase periodic multifunctional composites in three dimensions. One phase is stiff 
and insulating whereas the other is conductive and mechanically compliant. The 
optimization objective is to maximize a linear combination of the effective bulk 
modulus and conductivity of the composite. Composites with the Schwartz primitive 
and diamond minimal surfaces as the phase interface have been shown to have 
maximal bulk modulus and conductivity. Since these composites are not elastically 
isotropic their stiffness under uniaxial loading varies with the direction of the load. 
An isotropic composite is presented with similar conductivity which is at least 23% 
stiffer under uniaxial loading than the Schwartz structures when loaded uniaxially 
along their weakest direction. Other new near-optimal isotropic composites are 
presented, proving the capablities of the level-set method for microstructure design. 
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1 Introduction 



The design of multifunctional materials is a growing field. iTorquato et al. 
( I2OO2L I2OO3I ) considered the simultaneous transport of heat and electricity in 
three dimensional biphasic composites where one phase was more thermally 
conducting and less electrically conducting than the other. The structures were 
required to be isotropic by imposing threefold refiection symmetry. When both 
thermal and electrical conductivity were maximized the resultant composite 
resembled a Schwartz primitive (P) structureHThe optimality of the Schwartz 
P structure and also of the Schwartz D structure for maximizing both conduc- 
tivity properties was shown numerically via finite elemen t calculations and us- 
i ng rig orous cross-property bounds. Following this work. iTorquato and Donev 
( I2OO4I ) used finite element calculations to show the optimality of Schwartz P 
and D structures for maximum bulk modulus and conductivity for any ill- 
ordered phase properties (one phase has a larger conductivity but a smaller 
bulk modulus and shear modulus than the other phase). These publications 
demonstrated that structures derived from Schwartz P and D minimal surfaces 
are important multifunctional composites. 



Guest and PrevostI (120061 ) designed microstructures to maximize bulk modulus 
and fiuid permeability. They demanded cubic elastic symmetry and isotropic 
fiow synimetry and used the solid isotropic material with penalization (SIMP, 
Bends^d . Il989l : iBends^e and Kikuchil . 1 19881 ) implementation of topology opti- 
mization. Stokes' fiow conditions were assumed to calculate fiuid permeability 
and optimized structures were presented for different coefficients in the op- 
timization objective. For particular coefficients in the multiobjective design 
problem they obtained struc tures very similar to t he Sc hwartz P structure, 
consistent with the results of ITorquato and DonevI (120041 ). 



The recent paper of Ide Kruijf et al.l (120071 ) considered optimal structures with 
maximum stiffness and minimum resistance to heat dissipation and the de- 
sign of two dimensional composite materials with maximal effective thermal 
conductivity and bulk modulus. The two phases for the design problem were 
ill-ordered. The microstructures were required to be isotropic with respect to 
conductivity but only square-symmetric with respect to elasticity. The SIMP 
topology optimization algorithm was used. 



In this paper, two-phase isotropic three dimensional periodic composites are 
designed using topology optimization to have maximal bulk modulus and con- 
ductivity. The two isotropic phases for the microstructure design problem are 
ill-ordered infinite-contrast materials: one of the phases has finite stiffness 
but zero conductivity, whereas the other phase has finite conductivity and 



^ In this paper Schwartz P and diamond (D) structures refer to composites with 
the Schwartz P and D minimal surfaces as the phase interface. 
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zero stiffness. The result is competition between the phases. In particular, for 
the composite to have nonzero stiffness and nonzero conductivity both phases 
need to be connected. 



Our study is motivated by a desire to design maximally stiff, electrically con- 
ducting and isotropic cermets (composites of ceramic and metal), and opti- 
mally stiff, porous and isotropic bone implants. In the former situation, the 
ceramic has high stiffness but low conductivity, whereas the metal has compar- 
atively low stiffness and high electrical conductivity. For implants, the implant 
material (titanium, for example) has high stiffness but is impenetrable to bone 
in-growth, whereas the pore space has zero stiffness and allows bone in-growth. 
Here, the effective conductivity of the pore space is used to model the ease 
of bone in-growth into the implant. In reality, both scenarios will also involve 
manufacturing constraints, we deal with the simplest microstructure design 
problem and such constraints are left for future work. The calculation of elec- 
trical conductivity, thermal conductivity, dielectric constant, magnetic perme- 
ability and diffusion coefficient are all mathematically equivalent, thus our re- 
sults can also be applied to those cases. As well as having practical application, 
the objective function we consider is conveni ent from a theoretical viewpoin t 
due to the cross-property bounds derived by iGibiansky and Torquatd (119961 ) . 



The requirement of isotropy with respect to elasticity has not been seen in 
previous multifunctional composite design work, but can be an important cri- 
teria when the directions along which loads will be applied is unknown. In 
this case it is the strength of the weakest direction which is critical. In partic- 
ular, we highlight the anisotropic stiffness properties and weak directions of 
the Schwartz P and D structures to justify the isotropy constraint. As already 
noted, these structures were show n to be numerically optima l for maximum 
bulk modulus and conductivity by lTorquato and DonevI (120041 ). They possess 
only cubic symmetry so are not elastically isotropic and cannot be optimal in 
the present context. 



We employ the level-set implementation of topology optimization (I Allaire et al. 



2004l : IWang et al.U2003l ) to find optimized microstructures. Microstructure de- 
sign problems have not been attempted using this approach; we explore the 
capabilities of the method and develop a new way for imposing constraints. 



The remainder of the paper is as follows. Section [2] outlines the topology opti- 
mization problem, describes the relevant cross-property bounds and motivates 
the isotropy constraint. Section [3] describes our computational approach, pay- 
ing particular attention to how the isotropy constraint is imposed. Section H] 
presents our results. Section [5] highlights technical issues and summarizes the 
results with a "phase-diagram" which shows how competition between the 
multiple objectives effects the topology of the optimized structures. Conclud- 
ing remarks are given in Section [61 
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2 Problem outline and motivation 

2.1 Problem description 

The phase properties used throughout this paper are Ei = 1, z^i = 0.3, ai = 0, 
^2 = 0, = and ^2 = 1, where Ei is the Young's modulus of phase Ui is the 
Poisson's ratio of phase i and is the conductivity of phase i. Throughout this 
article the subscript 1 will refer to the stiff, insulating phase and the subscript 
2 will refer to the conductive, compliant phase. These are also referred to as 
the "stiff" and "conducting" phases respectively. 

The microstructure is represented by a unit cell cube which may be peri- 
odically extended along each coordinate direction. All computational results 
presented here used a unit cell represented by 40 x 40 x 40 voxels. We minimize 
the objective function 

J = -( — + — a ), (1) 

where ^* and a* correspond to the effective bulk modulus and conductivity of 
the material microstructure and cj^ and cUcr are weights which can be chosen 
to dictate the relative importance of the two objectives in our multiobjective 
design problem. Throughout this article cD^ and cDcr denote ^ and ^ respec- 
tively. 

The microstructure is required to be isotropic with respect to both stiffness 
and conductivity. Specifically, the effective elasticity tensor A^jj^i must be of 
the isotropic form 

^ij/cT = ^^^ij^ki + {5ik5ji + 5ii5jk) , (2) 

where /x* is the effective shear modulus of the composite and A* = /^* — 

is the other effective Lame constant of the composite. Similarly, the effective 

conductivity tensor K^ - of the composite must be of the isotropic form 

^r°=^*%- (3) 

The volume fraction of each phase is fixed. Throughout this paper Vi refers to 
the volume fraction of the stiff phase and the volume fraction of the conductive 
phase is V2 = 1 — 

Details regarding finding effective properties and imposing the required con- 
straints are left to Section [31 
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2.2 Cross-property bounds 



Cross-property bounds for conductivity and bulk modulus rigorously restrict 
the effective properties of a composite to be in some allowable region of the 
a* — /^* plane. They are used here to demonstrate the near-optimality of our 
optimized microstructrures. 

The tightest known cross-property bounds between bulk modulus and conduc- 
tivity for three dim ensional two-phase isotro p ic or cubic-symmetric compos- 
ites were derived by Gibiansky and Torquato ( 19961 ). For the particular phase 



properties used in this paper the set of all possible (a*, /^*) pairs is bounded 
by the two lines between the points 

{(0,0),(aH5,0)} and {(0, 0), (0, k^s)} (4) 



and by the hyperbola parameterized by: 



(j* = (l --f)aHS - 



7(1 - 7)^2 



HS 



Icfhs - (1 - Vi)a2 



7(1 -7)«:- 



HS 



(1 - 'y)KHS - ViHi ' 



(5) 



where 7 G [0,1]. (Jhs and Km r denote the upper Hashin-S htrikman bounds 
on the effectiv e conductivity (|Hashin and Shtrikmanl . 1 19621 ) and the effective 
bulk modulus (iHashin and Shtrikman . 19631 ). 



2.3 Schwartz primitive and diamond structures 



Composites with Schwartz P and D minimal surfaces as the phase interface 
must have ef fective properties whi c h sati sfy the cross-property bounds above 



for Vi = \. iTorquato and DonevI (120041 ) showed via numerical calculations 



that Schwartz P and D structures have properties which lie on a particular 
point on the cross-property bounds for any ill-ordered phase properties. For 
the particular phase properties considered in this paper, the Schwartz P and 
D structures have effective properties numerically equal to (a^, ^^) = (|, 



3' 63^ 



Fig. [T] shows the cross-property bounds, the optimal point (a^, /^^), and the 
calculated properties of 40 x 40 x 40 voxel approximations to the Schwartz P 
and D structures for our particular phase properties. As expected, our calcu- 
lated points do not exactly coincide with the theoretically optimal point due 
to the use of only 64, 000 voxels. 
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7 = 1 in Eq. 5 




a* 7 = in Eq. 5 

Fig. 1. Effective properties of the 64,000 voxel Schwartz P and D structures, the 
theoretical optimal point (cr^,A^^) and the relevant cross-property bounds. The 
effective properties must lie within the region bounded by hyperbola shown and the 
cr* = and = axes. 

The conductive phase of the approximate Schwartz P and D structures is 
shown in Figs. [2] and [3] respectively. These figures also show the directional 
dependence of the effective Young's modulus for the Schwartz P and D struc- 
tures. 

The effective Young's modulus £"* in a particular direction can be measured 
by loading the material uniaxially along that direction. £"* is then the value 
of the applied stress divided by the resulting strain as measured along the 
loaded direction. Once the effective elasticity tensor for the composite has been 
calculated its Young's modulus can be readily determined for any direction. 
We define E^-^ and E^^^ as the value of the effective Young's modulus along 
the directions for which it is smallest and largest respectively. 

Table [H tabulates the calculated properties of 64, 000 voxel approximations of 
the Schwartz P and D structures. 

Table 1 



Effective properties of the approximate Schwartz P and D structures pictured in 
Figs. El and Bl 



Structure 


a* 




A 


E* ■ 

mm 


^max 


Schwartz P 


0.3288 


0.1542 


0.1925 


0.1918 


0.3040 


Schwartz D 


0.3302 


0.1551 


0.1635 


0.1800 


0.2795 



Notes: A is the anisotropy of the structure, as defined in Section |31 
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Fig. 2. The conductive phase of the 64, 000 voxel approximation to a Schwartz P 
structure (left), and the directional dependence of its effective Young's modulus 
(right). 




Fig. 3. The conductive phase of the 64, 000 voxel approximation to a Schwartz 
D structure (left), and its calculated effective Young's modulus for all directions 
(right). 

2.4 Motivation for the isotropy constraint 

Figs. [2] and [3] clearly show the anisotropic stiffness properties of the Schwartz 
P and D structures. This is also demonstrated by the large difference between 
E^-^ and E^^^ for the two structures, as tabulated in Table [H In many engi- 
neering applications the presence of weak directions would not be favourable. 
We note that the Schwartz P and D structures have cubic symmetry and 
therefore are isotropic with respect to conductivity. Structures without cubic 
symmetry have directions with low conductivity, this would also not be desir- 
able in some cases. An isotropy constraint is an obvious suggestion to ensure 
that the Young's modulus and conductivity are the same in all directions. In 
such a case the picture corresponding to the right hand picture in Figs. [2] and 
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[3] would be a sphere with radius 

More precisely, the two following intuitive result holds. 



24.1 Result 

Given isotropic and anisotropic microstructures with the same average value 
of any directional property X*, the isotropic microstructure has the largest 

min' 



24.2 Proof 

The property X* is averaged over all directions for the anisotropic microstruc- 
ture to give X*: 

27r TT 

^* = J^j J X*{9,<P)sm{ci>)dcPd9, (6) 

0=0 (f)=0 



where 9 and (f) are the standard azimuthal and polar angles respectively. 
Since X* is the average of X* for all directions and the microstructure is not 
isotropic, there must exist a direction for which X* < X*. Thus X^^^^^^ < X*. 

Since the isotropic microstructure has this same X* and is isotropic, X*'^^^ = 
X* for all directions. Obviously X^^^ = X*. 

We obtain X^'^^'''' < X* = X^^^^ proving the result. 



3 Topology optimization algorithm 



3. 1 Calculation of effective properties 



We us e the standard elastic homogenization approach, for example see lGarboczi and Day 
(I1995I ) and the references therein, and choose to write it simply as 



^Ij — ^ijki^ki^ 



(7) 



where a*^- represents the effective stress tensor (not to be confused with a* 
which represents the effective conductivity) and eli represents the global strain 
tensor. The usual summation convention is used. The effective elasticity tensor 
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A'ljj^i can be calculated by applying six independent strain fields corresponding 
to the six independent components of the symmetric strain tensor. These 
strains are applied to the unit cell with periodic boundary conditions. The 
resulting stresses are found using the finite element method and are averaged 
over the voxels to give the components of the elasticity tensor. 

Similarly for the conductivity case, 

j; = k:jE;, (8) 

where J* is the effective current vector and £"* is the global applied electric 
field. The effective conductivity tensor K^j is calculated by applying three 
independent electric fields to the unit cell with periodic boundary conditions. 
The electric currents are computed by solving Laplace's equation with the 
finite element method and are averaged over the voxels in the finite element 
mesh to give the components of the conductivity tensor. 

We use the following definitions of the effective bulk modulus and shear mod- 
ulus: 



^* ^ 20 ^^^'^'^ ^ ^^'^^'^ ~ 30"^^'^'^'' ^^^^ 

which are invariants of the effective elasticity tensor and are certainly correct 
for the isotropic case in Eq. [21 

Similarly, the following is our definition of the effective conductivity 

^* = Ik;,. (11) 



3.2 Measuring anisotropy 



We measure elastic anisotropy as the "distance" between the calculated A'^-^ 
and the "nearest" isotropic elasticity tensor 

y^ijki ^ijki J y^ijki ^ijki ) . . 

ijkl ijkl 

with A'l'^^i given in Eq. [2] using effective properties from Eqs. [9] and [IHl The 
denominator can be simplified to 9/^*^ + 20/i*^. 



A 



mech 



\ 
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Similarly, we measure conductive anisotropy as 



A 



cond 



\ 



*,iso 



(13) 



where K^-^^^ is as in Eq.[3]with a* defined in Eq.dH The denominator simplifies 
to 3a*2. 

As an overall anisotropy measure these two individual anisotropies are summed 
to give 

3.3 Isotropy constraints 

The numerators of ^^g^/i ^cond sums of squares. To impose ^ = 
a constraint is introduced for each term in the sum. For example, in the 
conductivity case the six constraints are Cp^p G {1, . . . , 6}, where 

A/3a*Ci = X*i - a*, 
V3a*C3 = X33-a*, 
V3a*C5 = A/2Xi*3, 

73^*^6 = A/2i^^2*3- (15) 

Only five of these are independent but we work with all six to explicitly retain 
symmetry in the numerical calculations. The factors of \/2 are included in (^4, 
C5 and Cq on the right hand side and the factors of y^a* are included on the 
left hand side to give A^^^^ = C^. 

Similarly 21 constraints Cp^p G {7, . . . , 27} can be written down for the elas- 
ticity case, 19 of which are independent. They are also normalized so that 
summing and squaring them gives ^^gc/i- 

3.4 Level-set method for topology optimization 

Topolo gical optimization is imp lemented usi n g the level-set approach formal- 



ized by I Allaire et al.l ( 120041 ) and IWang et al.l (120031 ) and based on earlier work 
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by lOsher and Santosal (l2001r ) and ISethian and Wiegmannl (120001 ) . The reader 



is referred to I Allaire et al.l (120041 ) for standard details of the level-set ap- 



proach; only significant deviations from the standard method a re mentioned 



here. N ote that the topological derivative, as used for example bv I Allaire et al. 



(|2005h . is not used in the work presented here. We have found that, as noted by 
Allaire et al.l (120041 ). topological changes are readily achieved using the level- 
set method in three dimensions without the topological derivative, rendering 
it unnecessary. 



In our implementation the level-set function is initialized using an approxi- 
mation to the signed distance function, found by finding the nearest voxel of 
opposite phase with respect to the standard metric on M^. This initialization 
procedure is also used to reinitialize the level-set function at each iteration of 
the algorithm. 

The Hamilton-Jacobi evolution equation is solved via a standard upwind 
scheme and the time-step for the numerical evolution must be less than that 
given by the Cou r ant-F riedrichs-Lewy (CFL) condition to ensure numerical 
stability (jSethianl . 1 19991 ). We choose to use a time-step 1% of the CFL value 
and do many level-set evolutions per iteration of the optimization algorithm, 
typically between 25 and 100. The number of level-set evolutions per iteration 
is gradually reduced as the optimization converges to a local minimum. 



The shape derivatives of the effective elasticity and conductivity tensor com- 
pon ents and K^j ar e readily computed from well-known shape derivatives 



see 



Allaire et al.l . l2004l . and the references therein) by utilizing the superpo- 



sition principle for strains and stresses in linear elasticity and electric fields 
and currents in conductivity. From these it is straightforward to calculate the 
shape derivatives of the constraint£i| and of the effective properties /^* and a* 
to give the shape derivative of the objective function. These are of the form 



{0) = -J e.mvp pe{i,...,27}. 



(16) 



Here Q is the region occupied by the stiff phase. It is being deformed by the 
map Xi ^ Xi + so the left hand side of each equation above is the standard 
shape derivative. The boundary of Q is dQ^ with outward normal n^. The 
quantities v and Vp are called the shape sensitivities of J and Cp. 



^ The scaling factors of the constraints from the denominators of the anisotropy 
measures are not considered as variables but as constant scaling factors for the 
purposes of this calculation. 
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3.5 Imposing the isotropy constraints 



In the situation with no constraints, 9i is often chosen to be 9i = vui^ and 
the normal velocity V of the phase interface is chosen as V = OiUi = v. 
This implements a steepest descent type algorithm under evolution with the 
Hamilton-Jacobi equation. A modification of this for t he constrained case is 



outlined below. A more detailed presentation is given in lWilkins et al.l (120071 ). 



At any generic iteration of the algorithm, the constraints will be nonzero. To 
reduce both the objective function and the constraints we can require: 



J ^^n^'U = maximum possible, (17) 

on 

I e,n,Vp = rCp pG{l,...,27}, (18) 



on 



where r dictates the rate of exponential decay of the constraints. 

The shape sensitivities are functions defined on dQ and may be thought of as 
elements of a Hilbert space. In the following, the notation 1 1 • 1 1 and (•, •) is short 
hand for the norm and inner product on dQ: \ P = Jqq v^^ {vp^ Vq) = Jqq VpVq. 
In the finite element implementation Jqq VpVq is approximated by a sum over all 
the boundary voxels of the product VpVq. A voxel is determined as a boundary 
voxel if any of its 26 neighbours (including those across periodic boundaries) 
are of the opposite phase. 

Linearly dependent {vp} are removed from the set. We use the Gram-Schmidt 
procedure to build a mutually orthogonal set {vp : p G {1,...,24}} from 
{vp} which spans the constraint shape sensitivities. From this we can form the 
projection operator P which projects vectors in onto the space which will 
leave the constraints invariant. 

To implement Eq. [18] we take 



24 24 - 



where the ctp are real numbers which are chosen to solve the lower-diagonal 
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system 



rC24 



iPill 

(^^1,^2) 
pill 

(^1,^3) 
ll^ill 



ll^2|| 






1^3 1 



{vi,V24} {V2,V24,) {V3,V24,) 
\\vi\\ \\V2\\ \\V3\\ 








1^^241 



(20) 



In essence, this has decomposed OiUi into the sum of two parts. The first part 
is orthogonal to the shape sensitivities. The second part is a hnear combina- 
tion of the shape sensitivities. The former has been chosen to decrease the 
objective function as in Eq. f fTTI ). while the latter has been chosen to decrease 
the constraints as in Eq. f fTSl ). 

The parameter r is chosen to ensure 1 — > and > a'^-^. Chang- 

ing a'^in eflFects how strongly the algorithm projects onto the constraints. Typ- 
ically we choose o^^^^ = 0.1. 



3.6 Imposing the volume constraint 



The volume constraint is imposed via the Karush-Khun-Tucker technique 
( IKarushl . Il939l : iKuhn and Tucked . Il95ll ) . In our setting this technique involves 
calculating the expected volume change under evolution and if necessary cor- 
recting the shape derivative to ensure the volume change will not make the 
phase volume fractions deviate further from their required values. This type of 
technique h as been used with the l evel-set method for topology optimization 
previously ( jWang and Wangl . l2006l ) . 



3. 7 Velocity extension via smoothing 



The above method for the isotropy constraints calculates the normal velocity V 
at all boundary voxels. The normal velocity is set to zero for all non-boundary 
voxels. 

One issue with the level-set method of topology optimization is that the ve- 
locity must be extended away from the boundaries to give evolution of the 
structure. Frequently this is addressed by using an ersatz material approach, 
whereby material properties of the weak phase are set to some small nonzero 
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value instead of to zero, see for example I Allaire et al.l (120041 ). We choose to 
apply a smoothing to the velocities via the convolution 



where the subscripts refer to indices of the voxels. This convolution is applied 
repeatedly until Hamilton-Jacobi evolution will result in a geometric change 
of the structure. We note that this smoothing operation allows us to use only 
two phases during the optimization: no intermediate densities are required. 



4 Results 



Any initial unit cell is suitable to start the optimization process provided the 
initial periodic material has nonzero stiffness and nonzero conductivity in all 
directions. All of the optimized structures presented are optimized from the 
initial unit cell shown in Fig. [4](a). 

The chosen initial unit cell has two nice properties. First, the two phases 
have the same topology and geometry, meaning that there is no initial bias 
towards a particular property. Although the two phases of the initial unit 
cell look different in Fig. [4](a), a shift of either phase along each coordinate 
direction by half the unit cell edge length demonstrates that the two phases 
of the initial microstructure are geometrically identical. Second, the cell has 
simple cubic symmetry. Our computational approach retains this throughout 
the optimization and Acond = at each iteration. Isotropy with respect to 
elasticity is more difficult to achieve. 

The choice of internal parameters in the algorithm effects the outcome of the 
optimization process. The differences are often insignificant and experience 
with the algorithm allows the user to find a small set of internal parame- 
ters to use on each optimization problem. Results presented reflect the op- 
timized structure with the best objective function value from optimizations 
with various choices of internal parameters. It is not surprising that other in- 
ferior structures can be reached with different parameter choices: the level-set 
method converges to a local minimum and we do not know that a single global 
minimum exists. 
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(a) (b) (c) (d) 

Fig. 4. An example optimization history with (lJi^^CJct) = (1, |)- The four columns 
are (a) the initial unit cell, (b) iteration 20, (c) iteration 55 and (d) the optimized 
unit cell. In each case the top and bottom picture show the stiff phase and conductive 
phase respectively. 

4^1 Equal phase volume fractions 

First we consider the design of unit cells with an equal volume fraction for each 
phase, i.e., Vi = ^ = V2. Fig. [4] shows an initial structure, two intermediate 
structures and the optimized structure for (a;^,6Jcr) = (1, |)- Fig. [4] highlights 
that no intermediate densities are used at any time in the optimization. 

The material properties for the four microstructures shown in Fig. [Hare given 
in Table [21 Note that the volume constraint and the isotropy constraint are 
not satisfied during the optimization process. The errors on these constraints 
are small for the optimized structure (and for other optimized structures pre- 
sented). We have seen that small geometric changes can produce large per- 
centage changes in A when it is small, with corresponding very small changes 
in V2 and the objective function. Thus we consider a value of ^ < 0.005 
to be small enough for the isotropy constraint to be satisfied. Also, a value of 

^yreqmred _ yactual^ ^ q QQQ5 considered gOOd enough given that the unit 

cells are represented with only 64, 000 voxels without intermediate densities. 

Particularly of interest in Table [2] is how the minimum and maximum Young's 
moduli E^-^ and E^^^ change throughout the optimization. As the optimiza- 
tion progresses and the anisotropy A decreases, the difference between the 
maximum and minimum Young's moduli also decreases. For the optimized 
structure we see that the maximum and minimum Young's moduli are less 
than 1% apart, highlighting the effectiveness of the isotropy constraint. This 
is the case for all optimized structures presented and therefore £"* = ^^^^^^ is 
tabulated in the remainder of the paper. 



15 



Table 2 

Effective properties of the four unit cells from the optimization history shown in 



Fig.H 


Case 


J 


a* 


K* 


Vi 


A 


E* ■ 

min 


7?* 

^max 


(a) 


-0.3016 


0.3170 


0.1432 


0.5000 


0.2826 


0.1444 


0.2951 


(b) 


-0.3063 


0.3261 


0.1432 


0.5109 


0.1806 


0.1901 


0.2885 


(c) 


-0.3170 


0.2987 


0.1676 


0.5194 


0.0909 


0.2322 


0.2902 


(d) 


-0.3187 


0.3226 


0.1574 


0.5003 


0.0029 


0.2373 


0.2389 



Optimized unit cells for different {uj,^^ujcr) pairs are shown in Fig. El The ef- 
fective properties for these optimized microstructures are given in Table [31 
Fig. [6] shows the effective bulk modulus and conductivity for the optimized 
structures along with the cross-property bounds for the correct volume frac- 
tions and phase properties. The optimized structures are very close to the 
cross-property bounds and structures at several points along the bounds were 
readily obtained by changing the coefficients cj^ and uOcr in the objective func- 
tion. We note that Fig[: . [5]fa) rese i nbles the isotropic, maximum bulk modulus 
structure presented by ISigmundl ( 120001 ). 



Table 3 

Effective properties of the five optimized structures with Vi — ^ shown in Fig. [5l 



Case 




J 


a* 




Vi 


A 


E* 


9 


(a) 


(1,0) 


-0.2231 


0.0000 


0.2231 


0.5000 


0.0003 


0.3039 





(b) 




-0.2170 


0.1174 


0.2053 


0.4999 


0.0014 


0.2857 


7 


(c) 




-0.2201 


0.2709 


0.1749 


0.4999 


0.0013 


0.2547 


10 


(d) 




-0.3187 


0.3226 


0.1574 


0.5003 


0.0029 


0.2379 


10 


(e) 


(0,1) 


-0.3971 


0.3971 


0.0000 


0.4996 


0.0000 


0.0000 






Note: g is the genus per unit cell, see Section [5l 



4-2 30% stiff volume fraction 



To further explore the capabilities of the level-set method for microstructure 
design and the near-optimal structures for the stiffness-conductivity problem, 
a study similar to the above was performed with a required stiff phase volume 
fraction oiVi = Optimized structures for = ^ and different coefficients 
in the objective function are presented in Fig. [71 The effective properties for 
these structures are given in Table [Hand are summarized alongside the correct 
cross-property bounds in Fig. [8l 
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(e) 

Fig. 5. Optimized unit cells with Vi = ^ and different weighting schemes in the 
objective function: (a) (c^^,cj^) = (1,0), (b) (cj^,c^a) = (1, jq), (c) (co^.tJa) = (1, |), 
(d) (CJi^^uJct) = (l?^)? {^) i^H.^^a) = (O?!)- Ill each case the left and right picture 
show the stiff phase and conductive phase respectively. 

Table 4 



Effective properties of the optimized structures with = ^ shown in Fig. [71 



Case 




J 


a* 




Vi 


A 


E* 


9 


(a) 


(1,0) 


-0.1140 


0.0000 


0.1140 


0.3000 


0.0022 


0.1459 





(b) 


(b^) 


-0.1105 


0.1004 


0.1085 


0.3000 


0.0004 


0.1383 


7 


(c) 


(b^) 


-0.1108 


0.2129 


0.1001 


0.3004 


0.0009 


0.1309 


7 


(d) 


(b^) 


-0.1243 


0.5163 


0.0727 


0.3000 


0.0039 


0.1006 


10 


(e) 


(b^) 


-0.3399 


0.5457 


0.0671 


0.2996 


0.0011 


0.0987 


10 


(f) 


(0,1) 


-0.6069 


0.6069 


0.0000 


0.2998 


0.0000 


0.0000 






4^3 70% stiff volume fraction 



This section presents optimized structures with a required stiff phase volume 
fraction of Vi = ^. Optimized structures are displayed in Fig. [9l Effective 
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0.1 0.2 0.3 0.4 

a* 

Fig. 6. Effective properties for the optimized structures with Vi — ^ given in Fig. [5l 
The curve indicates the cross-property bounds for the problem. 

properties for the optimized structures are presented in Table [5] and summa- 
rized alongside the cross-property bounds in Fig. [101 

Table 5 



Effective properties of the optimized structures with Vi — shown in Fig. [9l 



Case 




J 


a* 




Vi 


A 


E* 


9 


(a) 


(1,0) 


-0.3861 


0.0000 


0.3861 


0.6998 


0.0001 


0.5234 





(b) 


(1,5) 


-0.3808 


0.0819 


0.3536 


0.6999 


0.0006 


0.4940 


3 


(c) 




-0.3938 


0.1140 


0.3368 


0.7003 


0.0004 


0.4762 


3 


(d) 


(1,1) 


-0.4710 


0.1597 


0.3113 


0.7001 


0.0012 


0.4547 


10 



5 Discussion 

The optimized structures with Vi = ^ display the same set of topologies as 
the optimized structures with Vi = \. The optimized structures for both cases 
have effective properties very close to the cross-property bounds. 

As in Fig. [H it is frequently the conductive phase which extends out to make 
new topologies in the optimization process. Accordingly, it was more diffi- 
cult to obtain optimized structures close to the cross-property bounds for 
the ^1 = ^ case. We see familiar topologies from the other volume fraction 
cases for structures (a) and (d) in Fig. [9l However, along the low conductivity 
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(e) (f) 

Fig. 7. Optimized unit cells with = ^ and different weighting schemes in the 
objective function: (a) (c^^,c^^) = (1,0), (b) (c^^,c^^) = (1, ^), (c) (cj^,cj^) = (1, ^), 
(d) (cJi^.cJa) = (1^1^), (e) {(^K^(^a) = (1,^), (f) = (0,1). In each case the 

left and right picture show the stiff phase and conductive phase respectively. 



0.1 


1 


1 1 
(b) 


1 1 1 1 
• 

(c) 


1 1 1 1 1 














1 


1 1 


1 1 1 1 


1 1 1 1 V /iX 



0.1 0.2 0.3 0.4 0.5 0.6 

a* 



Fig. 8. Effective properties for the optimized structures with Vi = ^ given in Fig. [71 
The curve indicates the cross-property bounds for the problem. 

section of the curve the optimized structures ((b) and (c) in Fig. [9j) have a 
conductive phase which is in two disconnected parts. One part of the con- 
ductive phase provides effective conductivity whereas the other is present to 
satisfy the isotropy constraint. 

For all three volume fraction cases we have been unable to obtain structures 
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(a) (b) 




(c) (d) 

Fig. 9. Optimized unit cells with Vi — jq and different weighting schemes in the 
objective function: (a) (c^^,cj^) = (1,0), (b) (cj^,cDa) = (1, |), (c) (0,^,0^) = (1, ^), 
(d) {Oi^^Ocr) = (1, 1). In each case the left and right picture show the stiff phase and 
conductive phase respectively. 



0.4 




Fig. 10. Effective properties for optimized structures with Vi = ^ given in Fig. [9l 
The curve indicates the cross-property bounds for the problem. 

in the high-conductivity part of the cross-property bounds curve which have 
nonzero stiffness. For example, see the gap between structures (d) and (e) in 
Fig.O We believe this problem stems from the isotropy constraint: for a small 
effective stiffness the algorithm simply disconnects the stiff phase to make it 
trivially isotropic. Aiming for this section of the curve cDa- > UJ^ in the objective 
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function, thus zero stiffness does not affect the objective function significantly. 
Structures in this region may be found by using an objective function which 
includes a penalty for disconnecting each of the phases, e.g., J = ^ + 
Alternatively, significant experimentation with o^^^^ may give structures with 
the desired properties. 

For the Vi = jq case we were also unable to obtain a structure close to 
i^HSi^)- The geometry expected for this structure given those obtained for 
the other volume fraction cases appears impossible at l/i = ^. Swapping 
the phases of structure (a) from Fig. [71 which has = ^ and properties 
close to (0,^/^5-), would give a structure very close to {(Jhs^^) on the Vi = 
cross-property bounds. However, it is the isotropy requirement for the 
elastic phase which drives the algorithm to the body-centered arrangement 
of Fig. [71(a). Given that the algorithm moves towards disconnecting the stiff 
phase, the optimization is not driven towards the correct structure. Instead it 
tends towards the simple cubic structure seen at other volume fractions, and 
when it is unable to disconnect the stiff phase cannot then move towards the 
body-centered geometry. 

In Fig. [TT] we hypothesize the existence of optimal microstructures with the 
topologies of the optimized microstructures we have presented. The microstruc- 
tures are grouped based on their top ology using the genus g which represents 



the number of handles of an object ( jHydd . 1 19891 ). 



6 Concluding remarks 



This paper presents a method for material design using the level-set method of 
topology optimization. The method is utilized to design isotropic periodic com- 
posite materials in three dimensions which are maximally stiff and conduct- 
ing from two ill-ordered, infinite contrast phases. The optimized microstruc- 
tures presented here have properties very close to the relevant conductivity 
bulk modulus cross-property bounds, proving the capabilities of the level-set 
method for microstructure design. 

Until now, the only known optimal single-scale microstructures for this prob- 
lem are cubic symmetric and therefore have weak directions. An important 
component of our method is the ability to improve upon this by imposing 
an isotropy constraint. The isotropy of the optimized microstructures makes 
them attractive for engineering applications: provided the shear modulus of 
the microstructures is high, the composites will have a high Young's modulus 
in all directions. The isotropy requirement means that almost all of our op- 
timized microstructures have not been presented previously. Our method can 
readily be applied to other constraints. 
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^ 9 




Fig. 11. The effective properties of optimal stiff and conducting microstructures for 
all volume fractions (above) and a parameterization of this surface using the stiff 
phase volume fraction Vi and the parameter 7 from Eq. [5] (below). In each case 
the black dots correspond to particular microstructures inferred from the optimized 
structures we have presented. The lower "phase-diagram" is divided into areas for 
which we hypothesize the existence of optimal microstructures with a particular 
topology. The three vertical dotted lines represent the three stiff phase volume 
fractions considered in this paper. Solid lines represent boundaries between regions, 
but should be considered as approximate. Dashed lines represent the extent of what 
we can infer from the optimized structures presented. Corresponding regions can be 
seen on the surface in the upper diagram. 
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For various volume fractions we hypothesize the existence of optimal single- 
scale microstructures with the topologies of those presented. We have not 
been able to find near-optimal structures for the high-conductivity nonzero 
stiffness section of the cross-property diagram. This is a result of the the 
isotropy constraint; instead of constraining the microstructures to be isotropic 
while being only weakly stiff the algorithm disconnects the stiff phase to make 
it trivially isotropic. 

The freedom to choose the initial unit cell has not been discussed. Different 
initial cells can give different prescribed symmetries and preliminary investi- 
gations demonstrate that different initial cells will almost certainly result in 
different optimized microstructures close to the cross-property bounds. Thus 
there are possibly other microstructures with similar properties to those pre- 
sented; this would provide the designer with freedom to choose the microstruc- 
ture to suit their purposes. However, given the closeness of our optimized 
structures to the cross-property bounds, other microstructures found could 
not outperform them by much more than a few percent. 
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